<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Chebychev design of an FIR filter given a desired H(w)</title>
<link rel="canonical" href="/Users/mcgrant/Projects/CVX/examples/filter_design/html/fir_chebychev_design.html">
<link rel="stylesheet" href="../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Chebychev design of an FIR filter given a desired H(w)</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% "Filter design" lecture notes (EE364) by S. Boyd</span>
<span class="comment">% (figures are generated)</span>
<span class="comment">%</span>
<span class="comment">% Designs an FIR filter given a desired frequency response H_des(w).</span>
<span class="comment">% The design is judged by the maximum absolute error (Chebychev norm).</span>
<span class="comment">% This is a convex problem (after sampling it can be formulated as an SOCP).</span>
<span class="comment">%</span>
<span class="comment">%   minimize   max |H(w) - H_des(w)|     for w in [0,pi]</span>
<span class="comment">%</span>
<span class="comment">% where H is the frequency response function and variable is h</span>
<span class="comment">% (the filter impulse response).</span>
<span class="comment">%</span>
<span class="comment">% Written for CVX by Almir Mutapcic 02/02/06</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% problem specs</span>
<span class="comment">%********************************************************************</span>
<span class="comment">% number of FIR coefficients (including the zeroth one)</span>
n = 20;

<span class="comment">% rule-of-thumb frequency discretization (Cheney's Approx. Theory book)</span>
m = 15*n;
w = linspace(0,pi,m)'; <span class="comment">% omega</span>

<span class="comment">%********************************************************************</span>
<span class="comment">% construct the desired filter</span>
<span class="comment">%********************************************************************</span>
<span class="comment">% fractional delay</span>
D = 8.25;            <span class="comment">% delay value</span>
Hdes = exp(-j*D*w);  <span class="comment">% desired frequency response</span>

<span class="comment">% Gaussian filter with linear phase (uncomment lines below for this design)</span>
<span class="comment">% var = 0.05;</span>
<span class="comment">% Hdes = 1/(sqrt(2*pi*var))*exp(-(w-pi/2).^2/(2*var));</span>
<span class="comment">% Hdes = Hdes.*exp(-j*n/2*w);</span>

<span class="comment">%*********************************************************************</span>
<span class="comment">% solve the minimax (Chebychev) design problem</span>
<span class="comment">%*********************************************************************</span>
<span class="comment">% A is the matrix used to compute the frequency response</span>
<span class="comment">% A(w,:) = [1 exp(-j*w) exp(-j*2*w) ... exp(-j*n*w)]</span>
A = exp( -j*kron(w,[0:n-1]) );

<span class="comment">% optimal Chebyshev filter formulation</span>
cvx_begin
  variable <span class="string">h(n,1)</span>
  minimize( max( abs( A*h - Hdes ) ) )
cvx_end

<span class="comment">% check if problem was successfully solved</span>
disp([<span class="string">'Problem is '</span> cvx_status])
<span class="keyword">if</span> ~strfind(cvx_status,<span class="string">'Solved'</span>)
  h = [];
<span class="keyword">end</span>

<span class="comment">%*********************************************************************</span>
<span class="comment">% plotting routines</span>
<span class="comment">%*********************************************************************</span>
<span class="comment">% plot the FIR impulse reponse</span>
figure(1)
stem([0:n-1],h)
xlabel(<span class="string">'n'</span>)
ylabel(<span class="string">'h(n)'</span>)

<span class="comment">% plot the frequency response</span>
H = [exp(-j*kron(w,[0:n-1]))]*h;
figure(2)
<span class="comment">% magnitude</span>
subplot(2,1,1);
plot(w,20*log10(abs(H)),w,20*log10(abs(Hdes)),<span class="string">'--'</span>)
xlabel(<span class="string">'w'</span>)
ylabel(<span class="string">'mag H in dB'</span>)
axis([0 pi -30 10])
legend(<span class="string">'optimized'</span>,<span class="string">'desired'</span>,<span class="string">'Location'</span>,<span class="string">'SouthEast'</span>)
<span class="comment">% phase</span>
subplot(2,1,2)
plot(w,angle(H))
axis([0,pi,-pi,pi])
xlabel(<span class="string">'w'</span>), ylabel(<span class="string">'phase H(w)'</span>)
</pre>
<a id="output"></a>
<pre class="codeoutput">
 
Calling SDPT3 4.0: 1199 variables, 321 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints = 321
 dim. of socp   var  = 899,   num. of socp blk  = 300
 dim. of linear var  = 300
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
    NT      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
 0|0.000|0.000|3.9e+02|2.4e+01|1.9e+05| 0.000000e+00  0.000000e+00| 0:0:00| chol  1  1 
 1|0.986|1.000|5.4e+00|3.0e-01|2.7e+03|-1.206201e-03 -3.513921e+01| 0:0:00| chol  1  1 
 2|0.995|1.000|2.6e-02|3.0e-02|4.4e+01|-8.587413e-04 -3.114402e+01| 0:0:00| chol  1  1 
 3|1.000|0.955|4.0e-08|9.4e-03|4.3e+00|-2.413584e-03 -4.280144e+00| 0:0:00| chol  1  1 
 4|0.955|0.914|3.1e-08|1.1e-03|1.3e+00|-2.495524e-02 -1.327064e+00| 0:0:00| chol  1  1 
 5|1.000|0.580|1.6e-10|4.7e-04|7.3e-01|-6.526725e-02 -7.919301e-01| 0:0:00| chol  1  1 
 6|0.501|1.000|7.3e-11|3.0e-06|5.7e-01|-3.874148e-01 -9.557116e-01| 0:0:00| chol  1  1 
 7|0.851|0.960|1.6e-11|4.1e-07|1.6e-01|-5.680266e-01 -7.295716e-01| 0:0:00| chol  1  1 
 8|0.854|1.000|4.2e-12|3.0e-08|2.8e-02|-6.822029e-01 -7.098235e-01| 0:0:00| chol  1  1 
 9|0.972|0.985|1.2e-12|3.4e-09|7.3e-04|-7.064224e-01 -7.071558e-01| 0:0:00| chol  1  1 
10|0.986|0.989|1.6e-14|3.4e-10|1.7e-05|-7.070901e-01 -7.071074e-01| 0:0:00| chol  1  1 
11|0.983|0.989|5.0e-16|4.7e-12|5.2e-07|-7.071063e-01 -7.071068e-01| 0:0:00| chol  1  1 
12|0.563|0.969|5.5e-12|1.1e-12|2.9e-07|-7.071065e-01 -7.071068e-01| 0:0:00| chol  1  1 
13|0.568|1.000|2.0e-11|1.1e-12|1.7e-07|-7.071066e-01 -7.071068e-01| 0:0:00| chol  1  2 
14|0.570|1.000|8.6e-12|1.6e-12|9.9e-08|-7.071067e-01 -7.071068e-01| 0:0:00| chol  1  2 
15|0.572|1.000|3.7e-12|1.7e-12|5.7e-08|-7.071067e-01 -7.071068e-01| 0:0:00| chol  1  2 
16|0.576|1.000|1.6e-12|1.0e-12|3.3e-08|-7.071067e-01 -7.071068e-01| 0:0:00|
  stop: max(relative gap, infeasibilities) &lt; 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 16
 primal objective value = -7.07106748e-01
 dual   objective value = -7.07106781e-01
 gap := trace(XZ)       = 3.29e-08
 relative gap           = 1.36e-08
 actual relative gap    = 1.36e-08
 rel. primal infeas (scaled problem)   = 1.57e-12
 rel. dual     "        "       "      = 1.00e-12
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 1.7e+00, 6.1e+00, 8.9e+00
 norm(A), norm(b), norm(C) = 8.4e+01, 2.0e+00, 1.8e+01
 Total CPU time (secs)  = 0.30  
 CPU time per iteration = 0.02  
 termination code       =  0
 DIMACS: 1.6e-12  0.0e+00  9.2e-12  0.0e+00  1.4e-08  1.4e-08
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.707107
 
Problem is Solved
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="fir_chebychev_design__01.png" alt=""> <img src="fir_chebychev_design__02.png" alt=""> 
</div>
</div>
</body>
</html>